
from scipy import argmin,nan,array,argmax,roll,zeros,shape,ones,diff,argsort,concatenate,shape,isnan,mean,unique,nanargmin,histogram,sort,float64,unique ,ceil,nansum
from pylab import find,subplot,clf,axes,rcParams,plot, xlabel,ylabel, imshow,gca,xlim,find,axes,twinx,legend,figure,bar,gcf,title,close
import commands


import tools.calc
#from IV import IV

from tools.plotroutines import *
class IETS():
  """
  IETS data analysis on data in IVs
  """
  def __init__(self,IVs):
    """
    init function
    """
    self.IVs=IVs
            
  def doReport(self,**kw):
    self.IVs.doReport(**kw)
    deeps,peaks = self.findpeaksmux(numdI=1,mux=0,nfitlength=35)
    close(7);figure(7,figsize=(14,12));
    for i in range(3):
      gaps = self.findgapspacing(peaks,filterpeak=5*(1+i))
      if len(gaps)>1:
        subplot(3,3,1+3*i);tmp=(gaps[:,3]+gaps[:,4])/gaps[:,4];hist(tmp[isnan(tmp)==False],bins=100);xlabel('Peak higths');ylabel('Counts')
        subplot(3,3,2+3*i);G = gaps[:,0]-gaps[:,2];
        if len(G[isnan(G)==False])>0:
          hist(G[isnan(G)==False],bins=100);xlabel('Gap size');ylabel('Counts')
        subplot(3,3,3+3*i);hist(gaps[:,0][isnan(gaps[:,0])==False],bins=50);hist(gaps[:,2][isnan(gaps[:,2])==False],bins=50);xlabel('peak positions mV');ylabel('Counts')
    mysavefig("results/%s/" % self.IVs.basename,"report_%i_gapanalysis" % (7))
    ndc_traces = find((self.I[1][0]<-nanmax(abs(self.I[1][0]))*1/3.).sum(1))
    if (len(ndc_traces)>1):
      clf();ax1=subplot(211);self.PlotIV(ndc_traces,dI=0,mux=0)
      subplot(212,sharex=ax1);self.PlotIV(ndc_traces,dI=1,mux=0);xlim(self.IVs.V.min()*1e3,self.IVs.V.max()*1e3)
      mysavefig(self.IVs.basename,'NDC_auto')
    

    
  def findpeaks(self,nfitlength=15,minsep=5e-3):
    (deeps1,peaks1) = self.findpeaksmux(nfitlength,mux=0,numdI=2)
    (deeps2,peaks2) = self.findpeaksmux(nfitlength,mux=1,numdI=1)
    self.deeps = array(list(deeps1)+list(deeps2))
    self.peaks = array(list(peaks1)+list(peaks2))
    self.reducepeaks(minsep)
    
  def reducepeaks(self,minsep=5e-3):
    RES = []
    for i in range(int(max(self.deeps[:,0]))):
      h=argsort(self.deeps[(self.deeps[:,0]==i)*(self.deeps[:,-1]==1)][:,1]) # arg sort peaks from trace 1 and mux=1
      if len(h)==0: continue
      deps = self.deeps[(self.deeps[:,0]==i)*(self.deeps[:,-1]==1)][h] # sort peaks
      k,l = find( diff(deps[deps[:,2]==0][:,1]) > minsep ), find( diff(deps[deps[:,2]==1][:,1]) > minsep ) # select those peaks wich have a 
      if len(k)>0 and len(l)>0:
        RES.append(concatenate((deps[deps[:,2]==0][k],deps[deps[:,2]==1][l]),0))
    self.deeps = concatenate(RES,0)
    RES = []
    for i in range(int(max(self.peaks[:,0]))):
      h=argsort(self.peaks[(self.peaks[:,0]==i)*(self.peaks[:,-1]==1)][:,1])
      if len(h)==0: continue
      deps = self.peaks[(self.peaks[:,0]==i)*(self.peaks[:,-1]==1)][h]
      k,l = find(diff(deps[deps[:,2]==0][:,1])>minsep),find(diff(deps[deps[:,2]==1][:,1])>minsep)
      if len(k)>0 and len(l)>0:
        RES.append(concatenate((deps[deps[:,2]==0][k],deps[deps[:,2]==1][l]),0))
    self.peaks = concatenate(RES,0) 
  def findpeaksmux(self,nfitlength=15,mux=1,numdI=1):
    V=self.IVs.V+self.IVs.VmuxOffset[mux]
    if isnan(self.IVs.Imean[mux][numdI]).all():
      self.IVs.__numDI__(nfitlength,mux=mux,dI=numdI)
    self.IVs.__numDI__(nfitlength,mux=mux,dI=numdI+1)
    I =  self.IVs.Imean[mux][numdI]
    dI = self.IVs.Imean[mux][numdI+1]
    t = zeros(shape(dI))
    dI[abs(dI)<1e-16]=0
    t[dI<0] = -1
    t[dI>0] = 1 
    Va=array([list((V[1:]+V[:-1])/2.)]*len(dI))
    tmp = array([range(len(V)-1)]*len(dI))
    trace = ones(shape(dI))
    trace[:,self.IVs.Ns[0]] = 0
    trace = trace[:,:-1]
    NB = array([range(len(dI))]*(len(V)-1)).T
    dt = diff(t,1)
    tmp1 = (dt>0)*(trace==0) + (dt<0)*(trace==1)
    peaks = array(zip(NB[tmp1],Va[tmp1],trace[tmp1],I[tmp1],[mux]*len(find(tmp1))))
    tmp1 = (dt<0)*(trace==0) + (dt>0)*(trace==1)
    deeps = array(zip(NB[tmp1],Va[tmp1],trace[tmp1],I[tmp1],[mux]*len(find(tmp1))))
    return deeps,peaks
  def analyseFoundPeaks(self,dI=1,mux=1):
    RES=[];RES2=[]
    if len(self.IETS_Peaks['peaks'])==0:
      return False
    Ns=unique(self.IETS_Peaks['peaks'][:,0])
    for n in Ns:
       # n = iv number
      I= self.IVs.Imean[mux][dI][n]
      Ip=tools.calc.correctbasline(self.IVs.V[self.IVs.V>0],self.IVs.Imean[mux][dI][n][self.IVs.V>0],pol=3,n=25,zerocor=0)/self.IVs.Imean[mux][dI-1][n][argmin(abs(self.IVs.V))]
      Im=tools.calc.correctbasline(self.IVs.V[self.IVs.V<0],-self.IVs.Imean[mux][dI][n][self.IVs.V<0],pol=3,n=25,zerocor=0)/self.IVs.Imean[mux][dI-1][n][argmin(abs(self.IVs.V))]
      for p in self.IETS_Peaks['peaks'][self.IETS_Peaks['peaks'][:,0]==n,:]:
        RES.append([mean(abs(p[1:])),self.IVs.Gfit[n]]+[Ip[nanargmin(abs(self.IVs.V[self.IVs.V>0]-p[i]))] for i in [1,3]]+[Im[nanargmin(abs(self.IVs.V[self.IVs.V<0]-p[i]))] for i in [2,4]])
        RES2.append([mean(abs(p[1:])),self.IVs.Gfit[n]]+[I[nanargmin(abs(self.IVs.V-p[i]))] for i in [1,3]]+[I[nanargmin(abs(self.IVs.V-p[i]))] for i in [2,4]])
    self.peakhights = array(RES)
    self.peakhights2 = array(RES2)
    return True
  def selectPeaks(self,tol=5e-3,maxpeaks=50,npeak=[4,4],usemuxstart=1,usebothMUX=False):
    """
    Select symetric peaks within tol of   4 present
    """
    try:
      find(self.IVs.Stable).any()
    except AttributeError:
      self.IVs.stdIVs()
    try:
      len(npeak)
    except TypeError:
      npeak = 2*[npeak]
    try:
      len(tol)
    except TypeError:
      tol = 2*[tol]      
    if not vars(self).has_key('peaks'):
      self.findpeaks()
    def getVtmps(N,mux,RES):
      return (self.peaks[(self.peaks[:,0]==N)*(self.peaks[:,2]==1)*(self.peaks[:,-1]==mux),1],
              self.peaks[(self.peaks[:,0]==N)*(self.peaks[:,2]==0)*(self.peaks[:,-1]==mux),1],
              self.deeps[(self.deeps[:,0]==N)*(self.deeps[:,2]==0)*(self.deeps[:,-1]==mux),1],
              self.deeps[(self.deeps[:,0]==N)*(self.deeps[:,2]==1)*(self.deeps[:,-1]==mux),1])
    
    def findsymetricpeaks(N,mux,RES,Vmean=None):
      if not (N==find(self.IVs.Stable)).any() or (sum((self.deeps[:,0]==N)*(self.deeps[:,-1]==mux))>maxpeaks) or  (sum(self.deeps[:,0]==N)==0) or  (sum(self.peaks[:,0]==N)==0):
        return RES
      else:
        tmpV1,tmpV,tmpV_deeps,tmpV1_deeps = getVtmps(N,mux,RES)
        if len(tmpV1)==0 or len(tmpV)==0 or len(tmpV_deeps)==0 or len(tmpV1_deeps)==0:
          return RES
        s=1
        for s in [1,-1]:
         for v in tmpV[tmpV*s>0] if Vmean is None else Vmean[Vmean*s>0]:
          closestsV = array([v,tmpV_deeps[argmin(abs(v+tmpV_deeps))],tmpV1[argmin(abs(v-tmpV1))], tmpV1_deeps[argmin(abs(v+tmpV1_deeps))]])
          meanV = sum(abs(closestsV))/4.
          if (sum((abs(closestsV)-meanV) < tol[mux]) >= npeak[mux]):
            (RES['peaks'] if s==1 else RES['dips']).append([N, tmpV[argmin(abs(meanV*s-tmpV))], tmpV_deeps[argmin(abs(meanV*s+tmpV_deeps))],tmpV1[argmin(abs(meanV*s-tmpV1))],tmpV1_deeps[argmin(abs(meanV*s+tmpV1_deeps))]])
        return RES
    
    RES={'peaks':[],'dips':[]}
    self.IETS_Peaks = []
    for N in range(len(self.IVs.Imean[0][0])):
      findsymetricpeaks(N,usemuxstart,RES)
    if len(RES['peaks'])==0 or len(RES['dips'])==0:
      print "no peaks found!!"
      self.IETS_PeaksTrace = {'peaks':[],'dips':[]}
      self.IETS_Peaks = {'peaks':[],'dips':[]}
      self.IETS_Peaks2=self.IETS_Peaks
      return None
    self.IETS_Peaks2 = {'peaks':tools.calc.deldup(array(RES['peaks'])),'dips':tools.calc.deldup(array(RES['dips'])) }
    RES2={'peaks':[],'dips':[]}
    if usebothMUX:
      for N in range(len(self.IVs.Imean[0][0])):
        Vmeans = list(abs(self.IETS_Peaks2['peaks'][:,1:]).mean(1))+list(- abs(self.IETS_Peaks2['dips'][:,1:]).mean(1))
        findsymetricpeaks(N,usemuxstart-1,RES2,Vmean=Vmeans)
      self.IETS_Peaks = {'peaks':tools.calc.deldup(array(RES2['peaks'])),'dips':tools.calc.deldup(array(RES2['dips'])) }
    else:
      self.IETS_Peaks = self.IETS_Peaks2
    self.IETS_PeaksTrace = {}
    for p in ['peaks','dips']:
      self.IETS_PeaksTrace[p] = array(map(lambda x: self.IVs.SavedTraces[find(self.IVs.SavedTraces[:,3]==x),4],self.IETS_Peaks[p][:,0])).ravel() if len(self.IETS_Peaks[p])>0 else []

  def __plotIETS_peak_hist__(self,trace=0,abstrue=False,peaktype='peaks'):
    Peaks = self.IETS_Peaks[peaktype][self.IETS_PeaksTrace[peaktype]==trace,1:].ravel()
    Peaks = Peaks[isnan(Peaks)==False]*1e3
    if abstrue:
      Peaks = abs(Peaks)
    return histogram(Peaks,bins=self.Vbin)
    
  def plotIETS_peak_hist(self,vbinnr=1,smooth=5,pertrace=True,datasets=None,minlentrace=10,figsize=(8, 6),**kw):
    if not vars(self).has_key('IETS_Peaks'):
      self.selectPeaks()
    V=self.IVs.V+self.IVs.VmuxOffset[1]
    Vbin = array(sort(dict(zip(V,V)).keys()))*1e3
    dV=Vbin[1]-Vbin[0]
    self.Vbin = array([Vbin[0]-dV/2.]+list((Vbin[1:]+Vbin[:-1])/2.)+[Vbin[-1]+dV/2.])
    self.Vbin = self.Vbin[::vbinnr]
    HISTs = {'peaks':[],'dips':[]}
    HISTs2 = {'peaks':[],'dips':[]}
    FigRes={};self.HISTsTRACES={};TRACES={};
    for peaktype in ['peaks','dips']:
      TRACES[peaktype] = unique(self.IETS_PeaksTrace[peaktype]) if datasets is None else datasets
      TRACES[peaktype] = filter(lambda i: shape(array(filter(lambda x: x[4]==i and not isnan(i),self.IVs.SavedTraces)))[0]>minlentrace,TRACES[peaktype])
      lTRACES = len(TRACES[peaktype])
      self.HISTsTRACES[peaktype]=TRACES[peaktype]
      FigRes[peaktype]={}
      for i in TRACES[peaktype]:
        HISTs[peaktype].append(self.__plotIETS_peak_hist__(trace=i,peaktype=peaktype)[0])
      HISTs[peaktype] = array(HISTs[peaktype])
      FigRes[peaktype][0] = figure(0+(0 if peaktype=='peaks'else 4),figsize=figsize);clf();
      if type(HISTs[peaktype].sum(0))==float64:
        print "HISTs (%s) is not type float" % peaktype
      elif type(self.Vbin[:-1])==float64:
        print "self.Vbin is not type float"
      elif len(self.Vbin[:-1])==len(HISTs[peaktype].sum(0)):
        Vb=self.Vbin[:-1]
        his=HISTs[peaktype].sum(0)
        bar(Vb,his,color='k')
        tmp = HISTs[peaktype].sum(0).cumsum();   
        tmp = (tmp[smooth:]-tmp[:-smooth])/float(smooth);    
        X = (self.Vbin[:-1][:-smooth]+self.Vbin[1:][smooth:])/2
        plot(X,tmp,'r:',lw=2)
        xlim(min(self.Vbin),max(self.Vbin))
        ylim(0,ceil(max(his[Vb>10])/5.)*5)
      else:
        print "len of Vbin = %i,HISTs=%i  " % (len(self.Vbin[:-1]),len(HISTs[peaktype].sum(0)))
      xlabel("Phonon energy / meV");ylabel('Counts')
      if pertrace:
        FigRes[peaktype][2] = figure(2+(0 if peaktype=='peaks'else 4))
        X = (self.Vbin[:-1]+self.Vbin[1:])/2
        k=ceil(lTRACES**0.5)
        l=ceil(lTRACES/k)
        Vb=self.Vbin[:-1]
        for i,his in enumerate(HISTs[peaktype]):
          subplot(l,k,i+1)
          bar(Vb,his,color='k')
          xlabel("Phonon energy / meV");ylabel('Counts');title(TRACES[peaktype][i])
          xlim(min(self.Vbin),max(self.Vbin))
          ylim(0,ceil(max(his[Vb>10])/5.)*5)
      for i in TRACES[peaktype]:
        HISTs2[peaktype].append(self.__plotIETS_peak_hist__(trace=i,abstrue=True,peaktype=peaktype)[0])
      HISTs2[peaktype] = array(HISTs2[peaktype])
      FigRes[peaktype][1] = figure(1+(0 if peaktype=='peaks'else 4), figsize=figsize);clf();
      #print shape(HISTs2[peaktype])
      if type(HISTs2[peaktype].sum(0))==float64:
        print "HISTs2 (%s) is not type float" % peaktype
      elif type(self.Vbin[:-1])==float64:
        print "self.Vbin is not type float"
      elif len(self.Vbin[:-1])==len(HISTs2[peaktype].sum(0)):
        Vb=self.Vbin[:-1]
        his=HISTs2[peaktype].sum(0)
        bar(self.Vbin[:-1],HISTs2[peaktype].sum(0),color='k')
        tmp = HISTs2[peaktype].sum(0).cumsum();   
        tmp = (tmp[smooth:]-tmp[:-smooth])/float(smooth);    
        X = (self.Vbin[:-1][:-smooth]+self.Vbin[1:][smooth:])/2
        plot(X,tmp,'r:',lw=2)
        xlabel("Phonon energy / meV");ylabel('Counts')
        xlim(0,max(self.Vbin))
        ylim(0,ceil(max(his[Vb>10])/5.)*5)
      else:
        print "len of Vbin = %i,HISTs2=%i  " % (len(self.Vbin[:-1]),len(HISTs2[peaktype].sum(0)))
      if pertrace:
        FigRes[peaktype][3] = figure(3+(0 if peaktype=='peaks'else 4),figsize=figsize)
        X = (self.Vbin[:-1]+self.Vbin[1:])/2
        Vb=self.Vbin[:-1]
        for i,his in enumerate(HISTs2[peaktype]):
          subplot(l,k,i+1)
          bar(Vb,his,color='k')
          xlabel("Phonon energy / meV");ylabel('Counts');title(TRACES[peaktype][i])
          xlim(0,max(self.Vbin))
          ylim(0,ceil(max(his[Vb>10])/5.)*5)
      if len(self.IETS_Peaks[peaktype])>0:
        print "Offset i V values of %s:%.2e" % (peaktype,(nansum(self.IETS_Peaks[peaktype][:,3])+nansum(self.IETS_Peaks[peaktype][:,4]))/float(len(self.IETS_Peaks[peaktype][:,4])))
      #print shape(self.IETS_Peaks[peaktype])
    self.HISTs=HISTs
    self.HISTs2 = HISTs2
    return FigRes
    
  def plotPeakIM(self,mux=1,dI=1,traces=None,normalize=False,**kw):
    self.IVs.PlotIM(mux=mux,dI=dI,normalize=normalize,reorderV=True,traces=traces,**kw)
    if len(self.IETS_Peaks)==0:
      return 0
    xt = xlim()
    yt = ylim()
    ax1=gcf().axes[0]
    if kw.get('single',False):
      ax4=axes([0.1,0.33,0.8,0.6],frameon=False,sharex=ax1);
    else:
      ax4=axes([0.1,0.33,0.8,0.3],frameon=False,sharex=ax1); 
      ax5=axes([0.1,0.63,0.8,0.3],frameon=False,sharex=ax1)
      if len(self.IETS_Peaks['peaks'])>0: ax5.plot(self.IETS_Peaks['peaks'][:,0]+.5,self.IETS_Peaks['peaks'][:,3:5],'kx')
      if len(self.IETS_Peaks['dips'])>0: ax5.plot(self.IETS_Peaks['dips'][:,0]+.5,self.IETS_Peaks['dips'][:,3:5],'k+')
      ax5.set_axis_off()
      ax5.set_ylim(self.IVs.V.min()+self.IVs.VmuxOffset[mux],self.IVs.V.max()+self.IVs.VmuxOffset[mux])
    if len(self.IETS_Peaks['peaks'])>0: ax4.plot(self.IETS_Peaks['peaks'][:,0]+.5,self.IETS_Peaks['peaks'][:,1:3],'kx')
    if len(self.IETS_Peaks['dips'])>0: ax4.plot(self.IETS_Peaks['dips'][:,0]+.5,self.IETS_Peaks['dips'][:,1:3],'k+')
    ax4.set_axis_off();ax4.set_ylim(self.IVs.V.min()+self.IVs.VmuxOffset[mux],self.IVs.V.max()+self.IVs.VmuxOffset[mux])
    xlim(xt);#ylim(yt)

  def printpeakstraces(self, DataSets=None, fileext='',minlentrace=1,closeIM=True,plotfunc = lambda : None,**kw):
    kw['single'] = kw.get('single',True)
    #print kw,self.HISTsTRACES
    self.plotPeakIM(normalize=True,**kw);mysavefig('results/%s/' % self.IVs.basename, 'IETS_IM_all_%s' % (fileext));close();
     #IETS_PeaksTrace = self.IETS_PeaksTrace if DataSets is None else self.IETS_PeaksTrace[find(map(lambda x: DataSets.count(x)==1,DataSetsAll))]
    HISTs2={};HISTs={};IETS_PeaksTrace={};DataSetsAll={};HISTsTRACES={}
    for peaktype in ['peaks','dips']:
      DataSetsAll[peaktype] = sort(dict(zip(self.IETS_PeaksTrace[peaktype],self.IETS_PeaksTrace[peaktype])).keys())
      IETS_PeaksTrace[peaktype] = self.IETS_PeaksTrace[peaktype] if DataSets is None else self.IETS_PeaksTrace[peaktype][find(map(lambda x: DataSets.count(x)==1,self.IETS_PeaksTrace[peaktype]))]
      HISTsTRACES[peaktype] = self.HISTsTRACES[peaktype] if DataSets is None else filter(lambda x: DataSets.count(x)==1,self.HISTsTRACES[peaktype])
      HISTs2[peaktype] = self.HISTs2[peaktype] if DataSets is None else self.HISTs2[peaktype][find(map(lambda x: DataSets.count(x)==1,self.HISTsTRACES[peaktype])),:]
      HISTs[peaktype] = self.HISTs[peaktype] if DataSets is None else self.HISTs[peaktype][find(map(lambda x: DataSets.count(x)==1,self.HISTsTRACES[peaktype])),:] 
      figure(figsize=[ 14,   6  ]);
      if len(HISTs2[peaktype])>0:
        bar(self.Vbin[:-1],HISTs2[peaktype].sum(0),color='k');xlim(0,None);xlabel("Phonon energy / meV");ylabel('Counts');xlim(0,max(self.Vbin));grid(True,lw=2);grid(True,which='minor');
        plotfunc()
      savefig('tmp_%s.png' % peaktype);close();
    plotIMs([(7/12.90,'results/%s/IETS_IM_all_%s.png' % (self.IVs.basename,fileext)),
               (6./14.,'tmp_peaks.png'),
               (6./14.,'tmp_dips.png')],width=20/2.)
    mysavefig('results/%s/' % self.IVs.basename, 'IETS_all_%s' % (fileext));close()
    for peaktype in ['peaks','dips']:
      figure(figsize= [ 14.2875,   7.2125])    
      if len(HISTs2[peaktype])>0:
        IMG = (HISTs2[peaktype][:,shape(HISTs2[peaktype])[1]/2.:].T.astype('f')/array(map(lambda x:sum(IETS_PeaksTrace[peaktype]==x),self.HISTsTRACES[peaktype]))).T
        IMG2 = (IMG.cumsum(1)[:,5:]-IMG.cumsum(1)[:,:-5])/5.
        imshow(IMG2,interpolation='nearest',extent=[0,max(self.Vbin),shape(HISTs2[peaktype])[0],0],aspect=max(self.Vbin)/shape(HISTs2[peaktype])[0]/1.5,cmap=my_cmap);colorbar();xlabel('meV');
        clim(0,0.1);ylabel('Trace')
        gca().xaxis.set_minor_locator(MultipleLocator(10));grid(True,lw=2);grid(True,which='minor')
      savefig('tmp_%s.png' % peaktype);close()
    ckw=copy(kw).tolist()
    dI=ckw.pop('dI',1)-1
    mux=ckw.pop('mux',1)
    self.PlotIM(normalize=True,dI=dI,mux=mux,**ckw);savefig('tmp2.png');
    if closeIM: close()
    plotIMs([  (7.2125/14.2875,'tmp_peaks.png'),
               (7.2125/14.2875,'tmp_dips.png'),
               (7/12.90,'tmp2.png')],width=20/2.)
    figtext(.5,1,self.IVs.basename,va='top',ha='center')
    mysavefig('results/%s/' % self.IVs.basename, 'IETS_IM2D_%s' % (fileext))
    if closeIM: close()
    files=['results/%s/IETS_IM2D_%s.pdf' % (self.IVs.basename,fileext),'results/%s/IETS_all_%s.pdf' % (self.IVs.basename,fileext)]
    DataSets = unique(list(self.HISTsTRACES['peaks'])+list(self.HISTsTRACES ['dips'])) if DataSets is None else DataSets
    for i in DataSets:
      L = array(filter(lambda x: x[4]==i,self.IVs.SavedTraces))[:,3]
      if len(L)<minlentrace:
        continue
      self.plotPeakIM(normalize=True,traces=list(L),**kw)
      mysavefig('results/%s/' % self.IVs.basename, 'IETS_IM_trace%i_%s'% (i,fileext))
      if closeIM: close()
      for peaktype in ['peaks','dips']:
        figure(figsize=[ 14,   6  ])
        if HISTsTRACES[peaktype].count(i)>0:
          j=HISTsTRACES[peaktype].index(i)
          bar(self.Vbin[:-1],HISTs2[peaktype][j],color='k');xlim(0,None);xlabel("Phonon energy / meV");ylabel('Counts');xlim(0,max(self.Vbin))
          gca().xaxis.set_minor_locator(MultipleLocator(10));grid(True,lw=2);grid(True,which='minor')
          plotfunc()
        mysavefig('results/%s/' % self.IVs.basename, 'IETS_Hist_trace%i_%s_%s'% (i,peaktype,fileext))
        close()
      figure(figsize=[ 14,   6  ]);title(i)
      plotIMs([(7/12.90,'results/%s/IETS_IM_trace%i_%s.png'% (self.IVs.basename,i,fileext)),
               (6./14.,'results/%s/IETS_Hist_trace%i_%s_%s.png'% (self.IVs.basename,i,'peaks',fileext)),
               (6./14.,'results/%s/IETS_Hist_trace%i_%s_%s.png'% (self.IVs.basename,i,'dips',fileext))],width=20/2.)
      mysavefig('results/%s/' % self.IVs.basename, 'IETS_trace%i_%s' % (i,fileext));close()
      files.append('results/%s/IETS_trace%i_%s.pdf' % (self.IVs.basename,i,fileext))
    commands.getstatusoutput('pdftk ' + (('%s '*len(files)) % tuple(files)) + (' cat output results/%sIETS_%s_collect.pdf' % (self.IVs.basename,fileext)))

  def printTraces(self,normalize=True, DataSets=None, fileext='',minlentrace=2,**kw):
    kw['single'] = kw.get('single',True)
    DataSetsAll = array(dict(zip(self.IVs.SavedTraces[:,4],self.IVs.SavedTraces[:,4])).keys())
    DataSetsAll = DataSetsAll[False== isnan(DataSetsAll)]
    self.IVs.PlotIM(normalize=normalize,**kw)
    autoclim(conf=.95)
    mysavefig('results/%s/' % self.IVs.basename, 'Traces_%s_IVs_all' % (fileext))
    DataSets = DataSetsAll if DataSets is None else DataSets
    files=['results/%s/Traces_%s_IVs_all.pdf' % (self.IVs.basename,fileext)]
    for j,i in enumerate(DataSets):
      L = array(filter(lambda x: x[4]==i,self.IVs.SavedTraces))[:,3]
      if len(L)<minlentrace:
        continue
      self.IVs.PlotIM(normalize=normalize,traces=list(L),**kw)
      autoclim(conf=.95)
      mysavefig('results/%s/'% self.IVs.basename, 'Plot_%s_IVs_tracse%i' % (fileext,i))
      close()
      files.append('results/%s/Plot_%s_IVs_tracse%i.pdf' % (self.IVs.basename,fileext,i))
    commands.getstatusoutput('pdftk ' + (('%s '*len(files)) % tuple(files)) + (' cat output results/%sIVs_%s_collect.pdf' % (self.IVs.basename,fileext)))
    return files